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ABSTRACT 

We investigate the impact of sinks of ionizing radiation on the reionization-era 21-cm signal, 
focusing on 1-point statistics. We consider sinks in both the intergalactic medium and inside 
galaxies. At a fixed filling factor of H II regions, sinks will have two main effects on the 21-cm 
morphology: (i) as inhomogeneous absorbers of ionizing photons they result in smaller and 
more widespread cosmic Hu patches; and (ii) as reservoirs of neutral gas they contribute a 
non-zero 21-cm signal in otherwise ionized regions. Both effects damp the contrast between 
neutral and ionized patches during reionization, making detection of the epoch of reionization 
with 21-cm interferometry more challenging. Here we systematically investigate these effects 
using the latest semi-numerical simulations. We find that sinks dramatically suppress the peak 
in the redshift evolution of the variance, corresponding to the midpoint of reionization. As 
previously predicted, skewness changes sign at midpoint, but the fluctuations in the residual 
HI suppress a late-time rise. Furthermore, large levels of residual HI dramatically alter the 
evolution of the variance, skewness and power spectrum from that seen at lower levels. In 
general, the evolution of the large-scale modes provides a better, cleaner, higher signal-to- 
noise probe of reionization. 

Key words: Key words: dark ages, reionization, first stars - intergalactic medium - methods: 
statistical - cosmology: theory. 


1 INTRODUCTION 

Reionization, the process in which the first galaxies ionized their 
surroundings and ultimately the entire Universe, is poorly con¬ 
strained at present. The 21-cm signal, as produced by a hyperfine 
transition in neutral hydrogen (HI), provides an excellent tracer of 
the neutral IGM. It is therefore the hope that future observations 
of the 21-cm brightness temperature using immense radio arrays, 
e.g. the LOw Frequency ARra>[^](LOFAR), the Murchison Wide- 
field Arra}0 (MWA), the Precision Array for Probing the Epoch 
of Reionizatioij^] (PAPER), the Hydrogen Epoch of Reionization 
Arra}j^](HERA) and the Square Kilometre Arra)[^](SKA), will dra¬ 
matically improve our understanding of reionization. 

To make optimal use of these observations it is vital that we 
build insight into the wide range of physical processes involved 
in reionization. Two of the most important factors to consider are 
the star formation history and the degree to which ionizing radi- 
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ation from galaxies escapes into the intergalactic medium (IGM). 
These two physical mechanisms have naturally been the focus in 


the modelling of reionization, with both se mi-numerical (e.g.j Zahn 


et al. 20071 Mesinger & Furlanetto 2007 

Choudhury et al. 

2009 

Mesinger et al. 20111 Santos et al. 2010, T 

lomas & Zaroubi 

2011 

and numerical simulations (e.g. Trac & Cen 2007 

Zahn et al. 

2007 

Baek et al. 20091 Trac & Gnedin 2011 [flliev et al. 

2014 ) being used 


to refine our understanding of both. These simulations describe 
how ionized regions form, grow and ultimately merge to ionize the 
entire IGM. However, it is not enough to understand the sources 
alone, we must also consider the sinks of ionizing radiation, i.e. the 
atoms that once ionized will recombine, wasting ionizing photons 
that would otherwise contribute to reionization. Sinks will have two 
main effects on the 21-cm signal: 

(i) as absorbers of ionizing photons they impact on both the tim¬ 
ing of reionization and the morphology of cosmic HII patches; 

(ii) as reservoirs of neutral gas they contribute a non-zero 21-cm 
signal in otherwise ionized regions. 

Sinks can reside in both the IGM and the interstellar 
medium (ISM), i.e. inside and outside of galaxies. In the ion¬ 
ized IGM at reionization redshifts, sinks mostly consist of 
diffuse structures and partially self-shielded gas clumps called 
Lyman-limit systems (LLSs, e.g. Furlan etto & Oh| ^2005 
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|McQuinn, Oh, & Faucher-Giguere|20lT] [Munoz et al.|2014|. 

LLSs are inhomogeneous in both space and time (Crociani et| 
[aT7||2011[ |Choudhury et al.||2009| »; as such, they can dramatically 
impact the duration and morphology of reionization (effect i above; 
|Sobacchi & Mesinger] |2014| . Although they are mostly ionized 
(e.g. jMcQuinn, Oh, & Faucher-Giguere||2011| |Rahmati et al.| 
|2013[ ), LLSs can also contribute to effect (ii). Indeed Sobacchi 
& Mesinger (2014) recently estimated that LLSs on average 
contribute a few percent neutral fraction inside cosmic H II patches 
during reionization. 


Sinks are also present inside galaxies where Hi can self 
shield in dense clumps. At lower redshifts (e.g. [Wolfe, Gawiser, &| 
|Proch aska 2005 ), sinks inside galaxies [mostly so-called Damped 
Lyman-alpha systems (DLAs); Prochaska et al. 2005, Noterdaeme 
|etal.|2 012 | provide a larger reservoir of HI compared to sinks in the 
IGM. Recent observations suggest that the contribution of DLAs to 
the cosmological mass density is roughly constant until it increases 
from z = 2.3 io z — 3.5 by almost a factor of 2 ( [Noterdaeme | 
|et al.||2012| ). There is currently no consensus in explaining these 
trends, making extrapolations to reionization redshifts uncertain. 
The abundance and ionization structure of such galactic sinks are 
more susceptible to the local environment, including stellar feed¬ 
back and baryonic cooling. Therefore it is difficult to quantify the 
contribution of galactic sinks to effect (ii) above. Using a ‘tuning- 
knob’ approach (similar to the one we use below), |Wyithe et al.| 
(20091 predict that by reducing the contrast between cosmic ion¬ 
ized and neutral patches (effect ii), galactic HI could indeed damp 
the 21-cm power spectrum by 10-20%. Galactic HI also contributes 
to effect (i) by reducing the efficiency with which ionizing radiation 
escapes into the IGM, usually quantified with the so-called ioniz¬ 
ing photon escape fraction, f esc . This clearly has an effect on the 
timing of reionization. However, reionization studies include the 
escape fraction and the associated uncertainties in the source term. 


In this work we examine the effect that both LLSs (IGM sinks) 
and differing levels of galactic Hi (galactic sinks) have on the 
statistics of the 21-cm brightness temperature, with emphasis on 
1-point statistics. We analyse the simulation presented in |Sobacchr| 
& Mesinger] < |2014| ), here after SM2014, to isolate the impact of 
IGM sinks on the 1-point statistics via effects (i) and (ii). We also 
test how the statistics of this simulation are altered in the presence 
of galactic sinks due to effect (ii). 


2 THE SIMULATIONS 


To investigate the impact of residual Hi on 21-cm statistics, we 
make use of the public simulation tool, 21CMFAS r l[^] The density 
and velocity fields in 21CMFAST are generated with standard per¬ 
turbation theory (Zel’dovich 1970). To generate ionization fields 
the excursion-set approach of Furla netto et al.| ( [2004| ) is applied to 
the evolved density fields: to determine whether a region is ionized, 
this algorithm compares at decreasing radii the number of ionizing 
photons to the number of baryons (plus recombinations). 

In this work, we study how sinks of ionizing photons impact 
the 21-cm signal, focusing on 1-point statistics. To do this, we run 
four 21CMFAST realisations of the ionization field during reioniza¬ 
tion, varying the impact of HI on both the large-scale reionization 
morphology and the level of residual H I inside cosmic H II regions 
(effects i and ii discussed in the introduction). We discuss the de¬ 
tails below, showing a summary of the four simulations used in 
Table [l] and Figure [l] All simulations are L = 300 Mpc on a side 
with a resolution of 400 3 . The density and velocity fields are the 
same in all of the runs. 


2.1 Homogeneous recombinations (21CMFASTvl) 

We call our reference simulation ‘NoLLS’. Sinks have a minimal 
impact in this simulation. It is generated with the current public re¬ 
lease, 21CMFASTvl. The resulting ionization fields are similar to 
ones obtained with radiative transfer simulations (e.g. |Zahn et ah] 
|2011| >. In 21CMFASTvl, a region is ionized if its (time-integrated) 
number of ionizing photons exceeds the number baryons plus a 
global average number of recombinations, n rec : 

C/coll(x, Z, 77, A7 m in) ^ 1 Tlrec 5 (1) 

where f co n(x, z, 77, M m - in ) is the fraction of matter, inside a region 
of scale 77, which resides in halos with mass greater than M m i n 
(taken here to correspond to a virial temperature of 10 4 K). The 
ionizing efficiency can be expanded as ( — /esc/*AT 7/ / b , where 
/esc is the fraction of ionizing photons that escape galaxies, /* is 
the fraction of galactic baryons inside stars, and N^/ b is the num¬ 
ber of ionizing photons produced per stellar baryon. This ionization 
criterion in eq. fT} is checked in an excursion-set fashion (F urlan-| 
etto et al.|2 004), starting from a maximum radius 77 ma x[^]down to 
the cell size. At the last smoothing scale, the unresolved, sub-cell 
Hu fraction is set to C(1 + rirec) _1 /coii(x, 2 ?, 77 ce ii, M min ). 

The ‘NoLLS’ run has a different reionization evolution and 
morphology from the other three runs [^] In order to facilitate direct 
morphological comparisons, we vary n rec in ‘NoLLS’ to match the 
redshift evolution of the filling factor of H II regions, Qmi from the 
other runs. 


The rest of this paper is structured as follows. In Section [2] 
we describe the simulations that we use to model reionization in¬ 
cluding sinks in both the IGM and in galaxies; in Section [3] we 
present our analysis methods and examine the 21-cm maps and 
their probability-density functions; we then consider the evolution 
of the 1-point statistics and our ability to constrain them with ra¬ 
dio telescopes present and future in Sections |3.4| and |3.5[ Sec¬ 
tion [4] considers how much the power spectrum might be altered 
in the presence of galactic sinks and finally in Section [5] we sum¬ 
marise our findings. We use the following cosmological parameters 
(Q m ,Q A ,Q b , h, cr 8 , n) = (0.28, 0.72, 0.046, 0.70, 0.82, 0.96), con¬ 
sistent with recent measurements by the Planck satellite ( [Planck| 
[Collaboration et al.|2014[ ). 


6 http://homepage.sns.it/mesinger/Sim.html 

7 T^max is a free parameter based on ionized photon mean free path at 
the redshifts of interest, see |Storrie-Lombardi et al.| |l994); Miralda-Escude 
(2003) ; [Choudhury et al. 1(2008) . In our ‘NoLLS’ run it is set to 30Mpc. 
However there is no maximum smoothing radius enforced by the simula¬ 
tions of SM2014 since the ionized photon mean free path is determined by 
the properties of ‘LLS’, which are computed self-consistently. 

8 The mean free path, 77 max , framework does not capture the sub-grid 
recombination physics of the SM2014 model, discussed in the next section. 
Nevertheless, these authors find that a choice of 77 max ~ 10 Mpc results 
in a reionization morphology and evolution roughly consistent with their 
sub-grid model. This ‘effective ionizing photon horizon’ is in fact much 
smaller than the physical mean free path (defined based on the instantaneous 
recombination rate and emissivity), and is empirically derived. 
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Table 1. Model summary. 


Model 

Properties 

( xhi)m in Hll 
regions at z = 

9.20, Qhh - 0.5 

NoLLS 

21CMFASTvl; sinks are in¬ 
cluded as a constant in the ionizing 
efficiency. Cosmic H n patches are 
fully ionized. 

0.0 

LLS 

Self-consistent treatment of local 

recombinations & UVB feedback 

of SM2014. Small absorbers act 
as sinks of ionizing radiation to 
slow the progress of reionization, 
reduce the size of ionized regions 
& suppress the 21 -cm temperature 
contrast between H n (over-dense) 
and H I (under-dense) regions. 

0.03 

LLS + 

Gal Hi 

‘LLS’ with the addition of a fixed 
percentage of remnant H I in galax¬ 
ies; its contribution to the 21-cm 
signal cosmic inside Hu regions 
further suppresses the contrast be¬ 
tween ionized and neutral patches. 

0.1/coli: 0.04 

0.5/coii: 0.07 

Mini + 

100% 

Gal Hi 

‘LLS’ with both sterile minihalos 
(cooled via H 2 ) as well as atomic 
cooling halos 100% neutral. De¬ 
scribes the most intense impact 
that remnant H 1 could have on the 

21-cm signal by suppressing the 
contrast between cosmic ionized 
and neutral patches. 

0.19 


integrating over the entire density distribution, the recombination 
rate per baryon in an ionized cell is: 


dn r < 


dt 




r+oo 

Jo 


Pv (A, z) AriHttB [1 - xm (A)] 2 dA 


(3) 

where Py (A, z) is the full (non-linear, sub-grid) distribution of 
overdensities: A = uh/tih (Miralda-Escude et al. 2000| >, and 
xhi (A) is equilibrium neutral fraction computed using the empir¬ 
ical self-shielding prescription of |Rahmati et al. | f20 13 1 > . Then the 
total, time-integrated number of recombinations per baryon, aver¬ 
aged over the local cosmic HII patch can be written as: 



dfl rec 

dt 



HII 


(4) 


where zm is the ionization redshift of a simulation cell. Photo¬ 
heating feedback from reionization (i.e. the depletion of the gas 
reservoir available for star-formation), is included in this simula¬ 
tion with a similar approach, in which the minimum halo mass ca¬ 
pable of retaining gas is computed using the local (cell’s) photo¬ 
ionization history fSobacchi & Mesinger||2013] >. Photo-heating 
feedback has a sub-dominant impact on reionization morphology, 
compared with inhomogeneous recombinations; however, since the 
sources and sinks are spatially correlated, the effects are additive in 
that they both prolong reionization history. 

‘LLS’ (corresponding to the ‘FULL’ simulation of SM2014) 
is our fiducial model for reionization history and morphology. The 
impact of inhomogeneous, sub-grid Hi on the large-scale reion¬ 
ization morphology can be estimated by comparing the ‘NoLLS’ 
and ‘LLS’ runs (see Figure [2]). Our ‘LLS’ run also includes resid¬ 
ual Hi inside the cosmic Hu patches, contributing on average a 
mass-averaged neutral fraction of a few percent inside HII regions 
(SM2014; see table 1). 



Figure 1. Evolution of mass-averaged ionized fraction for our 4 main mod¬ 
els (by construction, the volume-averaged evolution is the same in all mod¬ 
els). To track evolution of the morphology of Hu regions, we plot these 
four models as a function of the H II volume-filling factor hereafter. 


2.2 Inhomogeneous recombinations in the IGM 

As the next step, our 6 LLS’ run includes inhomogeneous recom¬ 
binations in the IGM using the model of SM2014 (their ‘FULL’ 
model). Such recombinations, driven by ~kpc scale gas clumps can 
slow the growth of large cosmic HII patches, primarily resulting in 
reduction of large-scale reionization structure. SM2014 implement 
these local recombinations by modifying the ionization criterion 
from eq. 0 to: 

C/coll[x, Z , R, AT m i n (x, z)] ^ 1 77, re c(x, z) . (2) 

In the model of SM2014, both the number of recombinations and 
minimum halo mass hosting galaxies depend on the local proper¬ 
ties of each ~ IMpc-scale reionization cell. Taking self-shielding 
into account via the parametrization of |Rahmati et al.| ( [2013| ), and 


2.3 Inhomogeneous recombinations in the IGM + galactic HI 


Our remaining runs use the same large-scale reionization morphol¬ 
ogy and redshift evolution as ‘LLS’ (see Figure [2}. However here 
we explicitly increase the level of residual Hi inside the ionized 
regions. Since the amount of HI inside galaxies is very difficult 
to predict, we adopt a ‘tuning-knob’ approach similar to |Wyithe| 


et al.|2 009 Specifically, we compute the averag^] collapsed frac¬ 


tion inside each cell, using the conditional (dependent on the cell’s 
density) mass fraction (e.g. |Bond et al. |1 99 1| [Lacey & Cole|1993|) 
normalized to the mean predicted by Sheth-Tormen (e.g. |jenkins 
et al.| < |2001| ); see also the discussion in |Barkana & Loeb||2004[ 


Mesinger et al. 2011). We then assign a fixed fraction, a, of the 


cell’s baryonic mass to be neutral, so that each cell’s galactic Hi 
number density is aA ce iinij/ C oii(M min , z, A ce ii, R ce ii). Keeping 


with the simplicity of the ‘tuning-knob’ approach, we use a uniform 
value of M m i n corresponding to a virial temperature of 10 4 K. We 
emphasise that this implementation considers only the additional 
signal that galactic sinks contribute; their impact on reionization 
timing is assumed to be accounted for in the ionizing efficiency of 
sources. 

We study different values of a; for the sake of brevity, we 
only present the simulations with 10% and 50% of the galactic 
mass as neutral hydrogen. To provide a sense of the range of val¬ 
ues that a might take we can extrapolate lower redshift constraints 


9 We ignore the Poisson noise associated with halo discreteness, since this 
affects scales too small to be observable with the SKA. 
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from DLAs. The contribution of galactic neutral hydrogen to the 
critical mass is found to be Q m ~ 10 -3 for z < 4 (Prochaska 


|et al.l[2005| |Noterdaeme et al.|2012) . This can be related to a as 
a — n H i/(fib/coii), which predicts a ~ 0.05 towards the end of 
reionization, i.e z ~ 6 , rising to a ~ 1 at z > 11 as / co ii decreases 
with increasing redshift. However at the redshifts for which Q m has 
been constrained, a < 0.03; as such, here after we will make the 
distinction between ‘reasonable levels of galactic sinks’, (a < 0 . 1 ), 
and ‘large levels’ (a > 0.5). 

It is important to note that the tail of the |Miralda-Escude et al.| 
( 2000 ) density distribution should, in principle, include gas inside 
galaxies. However, the distribution was calibrated to the lower over¬ 
density IGM, and has been shown to be a poor fit to the matter dis¬ 
tribution inside lower-redshift galaxies, even neglecting feedback 
(e.g.|Bolton & Becker| ( |2009| >; |McQuinn, Oh, & Fa ucher-Gi guere| 
( [201 if ). Moreover, the accuracy of the Miralda-Escude et al.| ( |2000| ) 
PDF is untested both at the high-redshifts (z ~ 10) of interest here, 
as well as in predicting the conditional (dependent on the local, 
large-scale overdensity) mass fraction. Hence in our approach we 
use the conditional, excursion-set approach, which is much bet¬ 
ter tested in predicting the collapse fraction. However, since the 
excursion-set galactic Hi is added on top of the SM2014 neutral 
fraction field, our approach is bound to double count some galac¬ 
tic H I, and the values of a should only be treated as approximate 
(lower limits). We find that due to the inside-out nature of the epoch 
of reionization (i.e. biased H II regions), the excursion-set approach 
generally predicts much higher collapsed fractions inside H II re¬ 
gions than MHR00 (see e.g. table 1); thus we do not expect the 
double-counting to have a large impact on our conclusions. 


Table 2. Instrumental specifications assumed for noise calculations. LO- 
FAR and SKA parameters are taken from |Mellema et al7|j2013) , HERA 
from Pober et al. ( 201 4}/ priv ate communications with A. Liu and MWA 
parameters from |Tingay et al.H2013) . 


Parameter 

MWA 

LOFAR 

HERA 

SKA 

Number of stations 
(iVstat) 

128 

48 

331 

450 

Effective area 

(Aeff/m 2 ) 

21.5 

804 

8.4e3/A/stat 

lO^/A/gtat 

Maximum baseline 

(firnax/m) 

2864 

3000 

360 

10 4 

Integration time 
(tint/hours) 

1000 

1000 

1000 

1000 

Bandwidth 
(Bf. MHz) 

6 

6 

6 

6 


3 BRIGHTNESS-TEMPERATURE MAPS & MOMENTS 
3.1 Computing ID statistics and noise 

Once the ionization boxes have been calculated as described, the 
brightness temperature (5Tb) can be computed according to: 


T —T 

5Tb — —^——^(1 — e~ 


1 + z 
T 

i 27. s 


T s 

1 + z 0.15 

10 Q m /i 2 


Xm(l + $) 


H(z)/(l + z) 


dv r / dr 

1/2 / n b fe 2 


(5) 


V 0.023 


inK, 


2.4 Inhomogeneous recombinations in the IGM + extreme 
galactic HI 

Finally, we consider an extreme model with the maximal amount 
of residual neutral hydrogen, ‘Mini+100% Gal Hi’. Following the 
same procedure as outlined in the previous section, we assign 100 % 
of the baryonic mass in 10 3 < T v i r halos as Hi. This model 
effectively assumes that minihalos with a virial temperature of 
10 3 < T v ir < 10 4 K, i.e. those that have cooled by molecular cool¬ 
ing, have been sterilized by a Lyman-Werner background so that 
they cannot form stars (e.g. Haiman, Abel, & Rees|2000]|Ricotti,| 
[Gnedin, & Shull|200T[ [Mesinger, Bryan, & Haiman|2006) . More¬ 

over, this model assumes that the time-scale of minihalo photo¬ 
evaporation is much longer than the duration of reionization (e.g. 
|Shapiro et al.||2004| |Iliev et aL||2005| |Ciardi et al.||2006| ). Our 
‘Mini+100% Gal Hi ’ run therefore serves to illustrate the most 
extreme impact that residual galactic Hi could have on 21 -cm ob¬ 
servations. 


where 5Tb is dependent on the overdensity 5 — p/~p — 1 , veloc¬ 
ity gradient dv r / dr and cosmology. The gas spin temperature T s 
measures the occupation levels of the two hyperfine energy lev¬ 
els involved in the Hi 21 -cm transition and determines the 21 -cm 
optical depth r UQ . It is assumed in this work that T s T 7 , i.e. 
T s is much higher than that of the cosmic microwave background 
T 7 . Whilst this is expected to be a reasonable approximation during 
reionization due to X-ray heating, the influence of spin temperature 
fluctuations may not be negligible in the early stages of reioniza¬ 
tion, i.e. Qhii <0.2 ([Pritchard & Furlanetto|2007|, Mesinger et al. 
|20l3l|Ghara et al.|2014>. 

In our analysis of the simulations described above we use 
the same approach as |Watkinson & Pritchard| ( |2014| ), hereafter 
WP2014. We measure the variance 62 and the skew S 3 of a simu¬ 
lated box according to 


& = jvr-El^-tfr b ] 2 , 

JVpix i =0 

1 Npix _ 

5 3 =* T7- J2 i 5Ti ~ <5Tb ] 3 ’ 

iVpiX z=0 


( 6 ) 


where 5Ti is the brightness temperature of the i th pixel, the 
average brightness temperature in a box is given by 5Tb — 
ATT* ^ 2 i=o X $Ti and /V p i x is the total number of cells in a box. 
We will consider two normalisations for the skew, the standard 
skewness 63/(62) 3//2 and the dimensional skewness 63/62 which 
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noLLS, z = 7.00, Qhi = 0.92 LLS, z = 7.00, Qhi = 0.92 a gal =0.1, 2 = 7.00, Qhi = 0.92 a gal =0.5, z = 7.00, Qhi = 0.92 minihalos, z = 7.00, Q m = 0.92 



noLLS, z = 9.20 ,Qhi = 0.50 LLS, z = 9.20,Q m = 0.50 a gal =0.1, z = 9.20 ,Qhi = 0.50 a gal =0.5, z = 9.20 ,Qhi = 0.50 minihalos, z = 9.20 ,Qhi = 0.50 




Figure 2. Brightness-temperature maps for 300 Mpc boxes normalised to z = 9; maps are resolved to pixels of (3 Mpc) 3 and presented with the pixel depth 
flattened into the page. Maps correspond to, from LLS (left), LLS + 10%GalH I (middle-left), LLS + 50%GalH I (middle-right) and minihalos + 100%GalH I 
(right). The top row is for z=7.0 when the H II volume filling factor is 0.92; the bottom row corresponds to z=9.2 at which time the H II volume filling factor 
is 0.50. Black regions are those with zero brightness temperature. 


was found to be a more natural normalisation for the skew of the 
neutral-fraction boxes of WP2014. 

We make order of magnitude calculations for the error induced 
by instrumental noise (using the instrument parameters in table [ 2 }, 
but do not consider foreground residuals. We assume the noise is 
independent on each pixel but identically distributed in a manner 
well described by Gaussian random noise with zero mean and stan¬ 
dard deviation cr^ o1se where 


2 

^noise 


2.9mK 


/ 10 5 m 2 \ 



2 


/l + z \ 4 ' 6 /71MHz lOOhours 

VlooV V ^ 


( 7 ) 


obeys the properties discussed above. We then use a test statistic 
for the m th moments, ATT* X^=o x ( x * — x z ) rn J to determine unbi¬ 
ased estimators for the moments of Equation [ 6 ] We then calculate 
the noise-induced variance on each unbiased estimator. These de¬ 
rived 1 -cr errors are then propagated onto the normalised skewness 
estimators 73 = S 3 /S 2 and 73 = *§ 3 /(*§ 2) 3 ^ 2 to give 


^§2 jy (^*^2(7 no i se + CLioise) : 


Vy, 


K 


^3 


where 


(S 2 

1 


l_i/. + (* 

2 v s 3 + 


(S 2) 3 


Vi 


s 3 


( S 2 ) 4 ' S2 


V, - 2 


S 3 

(s 2 y 


Cs 2 s 3 1 


( 8 ) 


-3^C& 


4 (S 2 ) 5 ' S2 SI S2S3 ' 


This expression is derived with the noise on the brightness tem¬ 
perature described by A T N = T sys /r]i \/Azz£int where the ar¬ 
ray filling factor is defined as rjf = ^tot/79max in which A to t 
is the total effective area of the array, D max = A/A 0\ the sys¬ 
tem temperature T sys is assumed to be saturated by the sky tem¬ 
perature at the frequencies of interest. Before measuring the mo¬ 
ments we smooth and resample the boxes to correspond roughly 
with the resolution expected from the various instruments; to do 
so, we use a Gaussian smoothing kernel with width R p i x = 6 Mpc 
for LOFAR/ MWA and i7 pix = 3 Mpc for HERA/ SKA (although 
we note that in this work, LOFAR/ MWA results are only pre¬ 
sented for larger smoothing scales). We also reduce the assumed 
frequency resolution of each telescope to correspond to R p i x , i.e. 
Av = H 0 v 0 VSUR pix /[cVTl + «)]• 

To estimate how the noise we’ve so far discussed will propa¬ 
gate onto the skewness and variance we assume that each pixel has 
a measured signal Xi — STi + m, where the noise on the pixel rii 


V§ = ~^—{3al oise K 4 + 12 S 2 cri oise + 5 < 7 ® oise ), (9) 

1 vpix 

and 

6 2 

C S 2 S 3 — ^^noise • (10) 

1 ' / P i X 

3.2 Brightness-temperature maps 

It is useful to take a look at the maps themselves before analysing 
their statistics. In the maps of Figure [2] ‘NoLLS’, ‘LLS’, ‘LLS + 
10%GalHl’, ‘LLS + 50%GalHl’ and ‘minihalos + 100%GalHl’ 
are represented from left to right; the top strip corresponds to z = 
7.0 and Q HU ~ 0.92 and the bottom strip corresponds to z = 9.20 
and Qhh 0.50. All maps have been smoothed and resampled to 
produce pixels of 3 Mpc on a side, the pixel depth has been flattened 
into the page. 
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As mentioned above, the impact of inhomogeneous IGM sinks 
can be seen by comparing the large-scale reionization morphology 
of the left panel to the other panels. As discussed in SM2014, re¬ 
combinations in unresolved structures result in a dramatic suppres¬ 
sion of large-scale HII regions. Comparing at a fixed Q hii, the dis¬ 
tribution of HII patch sizes is both narrower and shifted to smaller 
scales. The impact is dramatic enough that partially ionized cells, 
hosting unresolved Hu bubbles, become important. The ionized 
fraction of these cells is only <10%, but since they correspond to 
the densest (of those not yet reionized) ~Mpc-scale patches, this 
is sufficient to damp the signal from the neutral cosmic patches (as 
evidenced by the ‘yellow’ vs ‘orange’ IGM surrounding the HII 
bubbles in the ‘NoLLS’ and ‘LLS’ panels). This damps the con¬ 
trast between the fully ionized and fully neutral regions. We caution 
however that unresolved, sub-grid HII regions are taken into ac¬ 
count only in a simplistic fashion in our model, as described above; 
further work is required to accurately model such sub-grid physics. 

The impact of an increasing level of residual Hi inside the 
cosmic H II patches can be seen by comparing the insides of the 
‘black’ regions in the panels of the Figure [2] from left to right. The 
model of SM2014 (our ‘LLS’) only results in a few percent level 
of H I inside H II regions, and so on this colour scheme the H II 
regions are still black. However the extreme galactic HI models, 
‘LLS + 50%GalHl’ and ‘minihalos + 100%GalHT, show a no¬ 
table level of H I inside the cosmic H II regions, with small-scale 
structure tracing the distribution of galaxies. This further reduces 
the contrast between neutral and ionized regions, complicating an 
interferometric detection of reionization. In principle however, re¬ 
solving this small-scale Hi structure with high-resolution 21-cm 
interferometry could constrain models of galactic HI. 

3.3 Probability denstity functions 

These trends are quantified in Figure[3] where we present the prob¬ 
ability density functions (PDFs) from these maps, with lines cor¬ 
responding to: ‘NoLLS’ (pink dashed w/stars), ‘LLS’ (red solid), 
‘LLS + 50% galactic H T (black dotted w/triangles) and ‘minihalos 
+ 100% galactic Hi’ (green dashed w/circles). Plots correspond to 
Qhii ~ 0.3, 0.5, 0.7 (z =10.8, 9.2, 8), top to bottom. 

The presence of residual H I inside cosmic H II regions sup¬ 
presses zero-valued pixels in the PDFs. The bi-modal nature of the 
PDF is maintained for ‘LLS’ but the sharp peak at 5Tb = 0 mK 
is smeared into a wider peak centred on 5Tb ~ lmK because of a 
few percent level of residual Hi in the ionized IGM (SM2014). It 
is only with large levels of galactic sinks that the bi-modal nature 
is entirely suppressed, 

Also evident in these PDFs is the impact of sub-cell, unre¬ 
solved H II regions; these preferentially occur in the highest density 
pixels hosting newly forming galaxies, which is why there is a sup¬ 
pression of the high 5Tb tail in ‘LLS’ when compared to ‘noLLS’. 
This effect conspires with the loss of a sharp zero peak to further 
reduce the variance in the maps. 

3.4 Observing the brightness-temperature moments: 

Instrumental noise 

Figure [4] shows evolution of the variance (top), skewness (mid¬ 
dle) and dimensional skewness (bottom) as a function of the H II 
volume-filling factor. Lines correspond to ‘LLS’ (red solid), ‘LLS 
+ 10% galactic Hi’ (blue dot-dashed), ‘LLS + 50% galactic Hi’ 
(black dotted w/triangles), ‘minihalos + 100% galactic Hi’ (green 
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Figure 3. Brightness temperature PDFs at Qhii ~ 0.3, 0.5, 0.7 (z =10.8, 
9.2, 8), top to bottom. All maps are smoothed to a co-moving pixel size of 
3 Mpc. 

dashed w/circles) and ‘NoLLS’ (pink dotted w/ stars); we will use 
this key for the remains of the paper. 

We begin by reviewing what simulations that ignore recom¬ 
binations in unresolved systems predict for the evolution of the 
brightness-temperature moments. Studies of statistics from such 
models identified several important properties that would provide 
invaluable constraints on the timing of reionization if borne out in 
reality (e.g.|Furlanetto et al.|2004||Harker et al.|2009|| Watkinson &| 
Pritchard 2014). It is useful for the discussion that follows to refer 
to the pink dotted lines w/stars in Figure [4] The amplitude of the 
variance is maximised at the half-way mark of reionization as the 
weight of the non-zero distribution of the brightness-temperature’s 
PDF and its spike at zero are balanced (Figure[4|top). The late-time 
dominance of the zero-temperature spike induces a rapid increase 
in skewness at late times (Figure [4] middle). It could be argued that 
this effect is driven by the variance becoming very small as reion¬ 
ization draws to a close, however the dimensional skewness still 
exhibits a late-time signature but in the form of a turnover (Fig¬ 
ure [4] bottom). At the halfway point both skewness statistics pass 
from negative to positive. At early times they exhibit a minimum 
as the PDF transitions from being dominated by the density field 
to being dominated by the neutral fraction jLidz et al.|2008| ; the 
location of this minimum depends on the overlap of X-ray heating 
and reionization epochs (Mesinger et al. 2 013| ). 

Unresolved Hi sinks can dramatically change this picture. 
Already with ‘LLS’ the magnitude of the variance is suppressed 
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by up to a factor of two during the mid phases of reionization. 
Here the additional signal in H II regions, reduces the contrast be¬ 
tween ionized and neutral regions. More importantly, smaller H II 
regions mean that there is more sub-pixel, unresolved ionization 
structure in ‘LLS’ compared with ‘NoLLS’. This suppression of 
bi-modality (on resolvable scales) in the ionization structure is fur¬ 
ther exacerbated by the more disjoint Hu regions in ‘LLS’ which 
get ‘smeared-out’ when the maps are smoothed; i.e. the Hu re¬ 
gions in ‘LLS’ are less likely to be resolved by 21-cm interferom¬ 
eters. These effects dramatically suppress the small-scale contrast 
between ionized and neutral regions; the strength of the turnover 
associated with the halfway point of reionization is suppressed to 
the point that even SKA could have difficulty in constraining it. 
Fortunately, as we will see in the next section, further smoothing of 
the maps emphasises the turn-over feature, improving our ability 
to constrain this mid-phase signature. The turn-over in the variance 
occurs at Q H h ~ 0.4 in ‘LLS’ rather than nearly exactly halfway 
through reionization as predicted by ‘NoLLS’. This is because the 
ionization map is no longer well described as a two-phase field in 
which the number of pixels in each phase perfectly balance at the 
mid point of reionization. 

The inclusion of galactic HI reduces the variance even further 
(see also |Wyithe et al.|2009j ). Reasonable levels of galactic sinks do 
not dramatically alter the evolution of the variance (as could already 
be seen from the maps in Figure [2|. However, if the galactic Hi is 
tied to a fixed fraction of f co 11 as in our simple illustrative models, 
the growth of collapsed structure will boost the 21-cm variance at 
late times. In the extreme minihalo+100%Gal Hi model, structure 
growth totally dominates the signal throughout. SKA could easily 
detect such a feature. 

We see in the middle panel of Figure[4]that the transition from 
negative to positive skewness including sinks occurs at Q U u ~ 0.6 
(instead of Q m i ~ 0.5 as predicted by ‘NoLLS’). This transition is 
reduced to lower Q uu in the presence of extreme levels of galactic 
sinks. For reasonable levels of galactic sinks, this is predominantly 
due to effect (i), where small HII regions are unresolved. We even 
see this trend to a lesser degree when smoothing of ‘NoLLS’ on 
large scales smears-out the Hu regions (see Figure [6). For rea¬ 
sonable levels of galactic sinks, this signature and the turn-over in 
the variance bracket the mid-point and will constrain the timing of 
reionization more tightly than either statistic on its own. 

The late-time rapid increase predicted for the skewness is to¬ 
tally suppressed by the presence of IGM sinks; instead we see a turn 
over in the skewness of ‘LLS’ as the H II filling factor reaches about 
90%. This is because the normalising factor, i.e. the variance, is no 
longer tending towards zero and we instead see a transition from a 
skewness driven by the ionization field to one driven by the distri¬ 
bution of IGM sinks. This turn over is wiped out by the presence of 
even small amounts of galactic sinks. These are important points to 
be aware of as we begin our efforts to constrain reionization with 
21-cm observations. Detecting a turnover in the skewness would 
imply the presence of IGM sinks (as would reduced variance) and 
that reionization is in its final phases. The non detection of such a 
late-time signature in the skewness would not necessarily mean that 
we are not observing the end of reionization. Its absence combined 
with evolution in the variance with redshift would instead imply the 
presence of galactic sinks. 

SKA will be sensitive to the details of skewness we have de¬ 
scribed above, but only in the dimensional skewness, which ex¬ 
hibits more distinct features and better robustness to noise. 

The early-time minimum in the skewness statistics is seen to 
be robust across the models. This is understandable, since the lo- 


Redshift 

13.4 10.8 9.2 8.0 7.0 


120 


100 


b 80 


•e 

cd 

> 


40 


20 


8 

CO 


b 





- LLS 

. LLS+10% GalHl 

a-a LLS+50% GalHl 
• Mini+100% GalHl 

*•■■■★ NoLLS 
SKA error 


- LLS 

. LLS+10% GalHl 

a 'A LLS+50% GalHl 

• Mini+100% GalHl 

* * NoLLS 

SKA error 


. LLS+10% GalHl 

a 'A LLS+50% GalHl 
• Mini+100% GalHl 

*■■■■* NoLLS 
HI SKA error 


0.2 0.4 0.6 0.8 1.0 


HII filling fraction Qhii 


Figure 4. Variance (top), skewness (middle) and dimensional skewness 
(bottom) of brightness temperature measured from ‘observed’ maps with 
a co-moving pixel size of 3 Mpc as a function of volume-filling factor of 
ionized regions. Beige shadings depict SKA’s 1 -a instrumental errors for 
the ‘LLS’ simulation. Within the hatched region, spin-temperature fluctua¬ 
tions may not be negligible. 


cation of this feature is determined by the overlap of X-ray heat¬ 
ing and reionization ( [Mesinger et al.||2013| >, and here we assume 
Ts T 7 throughout. Nevertheless, extreme levels of residual HI 
will almost totally suppress this signature. At its default resolution 
SKA could struggle to detect this feature, depending on the tim¬ 
ings of reionization; however, its non-detection would still provide 
upper limits on this phase. 

Figure [5] shows the variance as a function of the percentage of 
galactic mass in sinks at the end of reionization, for this simulation 
z = 5.8. To reduce noise the maps are smoothed to a radius of 
Smooth = 10 Mpc and the beige shadings correspond to, from 
light to dark, the MWA, LOFAR and SKA 1-cr instrumental errors. 
Note that the signal when measured with SKA-like resolution, i.e. 
3 Mpc pixels, is greater by around a factor of ten. The blue dot- 
dashed line marks the level of residual hydrogen as constrained by 
DLAs at lower redshift, i.e. a = 0.03 at z ~ 4. 

There are two ways to view this post-reionization signal, one 
is that it is intrinsic noise to be overcome in our quest to constrain 
reionization, the other is that it provides a constraint on residual 
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Figure 5. Variance as a function of percentage of remnant galactic HI at 
z=5.80 (Qhii = 1) measured in ‘observed’ maps smoothed to a radius 
of ^smooth — 10 Mpc to reduce noise, the signal is around a factor of 10 
greater with the default pixel size of 3 Mpc. The vertical line corresponds to, 
a = 0.03, illustrating the level of residual hydrogen observed in galaxies 
at lower redshifts. MWA, LOFAR and SKA l-cr instrumental errors are 
depicted with beige shadings from lightest to darkest respectively. 


HI post-reionization. Clearly the amplitude of the variance at the 
end of reionization is very sensitive to the level of galactic sinks. 
This intrinsic noise will be detectable by SKA, which will therefore 
be well placed to constrain galactic sinks with this statistic. We do 
not explicitly plot HERA errors but they are indistinguishable from 
SKA, i.e. negligible, on this plot and so should perform equally 
well at this task. Before these telescopes come online, both LOFAR 
and MWA should be able to set some constraints on this quantity, 
although at the lower and more probable values of a, MWA’s sig¬ 
nal to noise will be too poor to do more than set upper limits on 
remnant galactic HI. The flip-side of this is that the intrinsic noise 
from galactic sinks should not be a limiting factor for these first 
generation instruments. 


3.5 Observing the brightness-temperature moments: 

Smoothing to reduce noise 

As in [Watkinson & Pritchard| ( [2014] ), we investigate the effects 
of smoothing the maps to beat down instrumental noise. Figure 
[6] (top) shows the evolution of the variance when smoothed to 
Smooth =10 Mpc as a function of Hll volume-filling factor. LO¬ 
FAR will be well placed to constrain the presence of IGM sinks 
through the reduction in the variance, could constrain/ exclude 
large levels of galactic sinks and will be able to offer timing con¬ 
straints from the mid-point turnover in the variance (unless the 
model is extreme and this signature does not exist). We see that 
smoothing recovers the late-time turnover in both skewness statis¬ 
tics for reasonable levels of galactic sinks. This is fortunate as the 
late-time maximum in the dimensional skewness will be detectable 
even by LOFAR for a wide range of galactic HI levels. The mid¬ 
point of reionization will be bounded by SKA and HERA measure¬ 
ments of the skewness, which will inform us when Q HU <0.7, and 
constraints from the variance of Q H ii ~ 0.4. The early-time min¬ 
imum should still be detectable in both skewness statistics despite 
suppression by the presence of IGM sinks. SKA is the only instru¬ 
ment that is likely to actually detect this feature; however LOFAR 
and HERA should be able to use its absence to set upper-limits on 
this phase. If residual hydrogen levels are extreme then this signa¬ 
ture, like every other, will be suppressed. 
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Figure 6. Variance (top), skewness (middle) and dimensional skewness 
(bottom) of brightness temperature as a function of the H II volume-filling 
factor measured in ‘observed’ maps smoothed to R sm 00 th =10 Mpc. Beige 
shadings depict l-cr instrumental errors for the ‘LLS’ simulation; tones cor¬ 
respond to MWA, LOFAR and SKA from lightest to darkest. Within the 
hatched region, spin-temperature fluctuations may not be negligible. 


4 BRIGHTNESS-TEMPERATURE POWER SPECTRUM 

The brightness temperature’s 1-point statistics form the focus of 
this paper, but the 21-cm power spectrum will be an incredibly im¬ 
portant statistic when it comes to constraining reionization. As such 
it is worth understanding the impact of galactic neutral hydrogen on 
this statistic. We concentrate on the spherically-averaged dimen¬ 
sionless power spectrum which we describe using A 2 Th (k,z) — 
k 3 J(27r 2 V)(\62i(k,z)\ 2 ) k in which S 2 i = 6T h (k,z)/6T h (z)-l, 
STb(z) is the redshift dependent average brightness temperature 
calculated from the simulation, V is the volume of the simulated 
box and the angle brackets denote an average over /c-space. To 
model instrumental errors on the spherically averaged power spec¬ 
trum, we take the approach outlined in the appendix of |McQuinn| 
|et al.| ( [2006] ) adopting a logarithmic bin width of e = 0.5. In calcu¬ 
lating the number density of baselines we assume a filled nucleus 
limited by the minimum separation of stations, followed by 

an r~ 2 drop off in station density. This assumes a smooth density of 
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k (Mpc -1 ) 

Figure 7. Dimensionless power spectrum of the brightness temperature as 
a function of fcMpc -1 for z = 10.8 (top), z — 9.2 (middle), and z — 8 
(bottom). 


stations, which for LOFAR is a particularly crude approximation. 
Whilst the exact configuration of SKA is yet to be decided, both 
LOFAR and MWA have optimized station positioning, as such our 
power spectrum errors are indicative only. 

The dimensionless power spectra for our simulations are pre¬ 
sented in Figure [7] for Qhii ~ 0.3, 0.5, 0.7 (z =10.8, 9.2, 8), 
from top to bottom respectively. As already discussed in SM2014, 
IGM sinks can suppress large-scale 21-cm power, resulting in much 
steeper power spectra throughout reionization. From Figure [7] we 
see that adding additional galactic HI does not have a large im¬ 
pact on the power spectrum, for reasonable values of a. More ex¬ 
treme levels of galactic Hi result in a further steepening of the 
power spectrum, since: (i) the contrast between ionized and neu¬ 
tral patches is reduced, and (ii) the distribution of galactic HI con¬ 
tributes additional power on small-scales. 

The evolution of the dimensional power spectrum 
5T b 2 (k,z) A 2 STb (k, z) as a function of the H II volume fill¬ 
ing factor is plotted in Figure [S] for k — IMpc -1 (top) and 
k = O.IMpc -1 (bottom). Unlike higher-order moments, the power 


Redshift 



HII filling fraction Qhii 

Figure 8. Dimensional power spectrum of the brightness temperature as 
a function of the Hu volume-filling factor for k = IMpc -1 (top) and 
k = O.IMpc -1 (bottom). Within the hatched region, spin-temperature 
fluctuations may not be negligible. 

spectrum at small k can cleanly pick-up large-scale structure. 
As such its timing signatures are more robust. For example, the 
mid-point turnover in the k = O.IMpc -1 dimensional power 
spectrum is robustly at Q H h ~ 0.5 for the more reasonable models, 
although it shifts to smaller Q mi with large levels of galactic sinks 
and is completely wiped out in the extreme galactic Hi model. 
We see that whilst reasonable quantities of galactic sinks have 
little qualitative impact, large quantities do impact on the nature 
of the power spectrum’s evolution with its evolution no longer 
well described by a simple inverted parabola on large scales. This 
results from additional small-scale power (evident in the plots of 
Figure [7] and the top plot of Figure [8j and from the fact that at late 
times there is still substantial residual HI that drives the power. 


5 CONCLUDING REMARKS 

In this paper we have examined the impact of sinks on the 21- 
cm statistics of reionization. We consider two classes of sink, IGM 
sinks and galactic sinks, which have two effects on the 21-cm sig¬ 
nal: 

(i) As absorbers of ionizing photons they impact on both the 
timing of reionization and the morphology of cosmic H II patches. 

(ii) As reservoirs of neutral gas they contribute a non-zero 21- 
cm signal in otherwise ionized cosmic patches. 

We use the simulations of |Sobacchi & Mesinger| ( |20T4| ) to inves¬ 
tigate the impact of IGM sinks on the 21-cm 1-point statistics due 
to effects (i) and (ii). In addition we study how the contribution of 
galactic sinks to (ii) affects the 21-cm statistics. We also consider 
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an extreme, illustrative model in which gas inside sterile minihaloes 
and atomically-cooled galaxies is fully neutral. 

It is likely that the impact of sub-grid IGM sinks on reion¬ 
ization morphology < |Sobacchi & Mesinger||2014] effect i), will 
suppress the expected ‘rise and fall’ reionization signal of the 
raw (unfiltered) 21-cm variance. Instead, the large-scale power (or 
smoothed variance) shows a more pronounced evolution during 
reionization, since it is less affected by unresolved H II structure. 
Thus first-generation instruments such as LOFAR should target 
large-scale modes for cleaner reionization probes. 

Similarly, the rapid increase of the skewness previously asso¬ 
ciated with the end of reionization is also suppressed by inhomo¬ 
geneous sinks. This is because the variance is non-zero even post¬ 
reionization (due to residual Hi; effect ii) and so we instead see a 
turnover as the skewness transitions from being dominated by the 
ionization field to the residual HI field. 

We find that more reasonable levels of galactic sinks, will not 
change the qualitative evolution of the moments during reioniza¬ 
tion. It is worth noting that the presence of even low levels (< 10%) 
of galactic sinks would result in the absence of any late-time skew¬ 
ness signature in SKA maps. However, the late-time turnover in the 
skewness should be recoverable by smoothing the maps to a lower 
resolution. 

We find that large levels (> 50%) of galactic HI have the 
potential to dramatically alter the qualitative evolution of moments. 
This provides an exciting opportunity for SKA (and even LOFAR) 
to constrain levels of galactic HI in the high-redshift Universe. 
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